This function takes the terms of a series in ARRAY of size
ARRAY_SIZE and computes the extrapolated limit of the series using
a Levin u-transform. The extrapolated sum is stored in SUM_ACCEL,
with an estimate of the relative error stored in PRECISION. The
actual term-by-term sum is returned in SUM_PLAIN. Additional
working space must be provided in the arrays Q_NUM, Q_DEN, DSUM of
size ARRAY_SIZE elements each and in DQ_NUM, DQ_DEN of size
ARRAY_SIZE**2. The algorithm estimates both truncation error and
round-off error to choose an optimal number of terms for the
extrapolation.
File: gsl-ref.info, Node: Example of accelerating a series, Next: Series Acceleration References, Prev: Acceleration functions, Up: Series Acceleration
Example of accelerating a series
================================
The following code calculates an estimate of \zeta(2) = \pi^2 / 6
using the series,
\zeta(2) = 1 + 1/2^2 + 1/3^2 + 1/4^2 + ...
After N terms the error in the sum is O(1/N), making direct summation
The ouput below shows that the Levin u-transform is able to obtain an
estimate of the sum to 1 part in 10^10 using the first eleven terms of
the series. The error estimate returned by the function is also
accurate, giving the correct number of significant digits.
bash$ ./a.out
term-by-term sum = 1.5961632439130233 using 20 terms
term-by-term sum = 1.5759958390005426 using 13 terms
exact value = 1.6449340668482264
accelerated sum = 1.6449340668166479 using 13 terms
estimated error = 0.0000000000508580
actual error = -0.0000000000315785
Note that a direct summation of this series would require 10^10 terms
to achieve the same precision as the accelerated sum does in 13 terms.
File: gsl-ref.info, Node: Series Acceleration References, Prev: Example of accelerating a series, Up: Series Acceleration
References and Further Reading
==============================
The algorithms used by these functions are described in the following
papers,
T. Fessler, W.F. Ford, D.A. Smith, HURRY: An acceleration
algorithm for scalar sequences and series `ACM Transactions on
Mathematical Software', 9(3):346-354, 1983. and Algorithm 602
9(3):355-357, 1983.
The theory of the u-transform was presented by Levin,
D. Levin, Development of Non-Linear Transformations for Improving
Convergence of Sequences, `Intern. J. Computer Math.' B3:371-388,
1973
File: gsl-ref.info, Node: Random Number Generation, Next: Random Number Distributions, Prev: Series Acceleration, Up: Top
Random Number Generation
************************
The GNU Scientific Library provides a large collection of random
number generators which can be accessed through a uniform interface.
Environment variables allow you to select different generators and
seeds at runtime, so that you can easily switch between generators
without needing to recompile your program. Each instance of a
generator keeps track of its own state, allowing the generators to be
used in multi-threaded programs. Additional functions are available
for transforming uniform random numbers into samples from continuous or
discrete probability distributions such as the Gaussian, log-normal or
Poisson distributions.
* Menu:
* General comments on random numbers::
* The Random Number Generator Interface::
* Random number generator initialization::
* Sampling from a random number generator::
* Auxiliary random number generator functions::
* Random number environment variables::
* Saving and restoring random number generator state::
* Random number generator algorithms::
* Unix random number generators::
* Numerical Recipes generators::
* Other random number generators::
* Random Number Generator Performance::
* Random Number References and Further Reading::
* Random Number Acknowledgements::
File: gsl-ref.info, Node: General comments on random numbers, Next: The Random Number Generator Interface, Up: Random Number Generation
General comments on random numbers
==================================
In 1988, Park and Miller wrote a paper entitled "Random number
generators: good ones are hard to find." [Commun. ACM, 31, 1192-1201].
Fortunately, some excellent random number generators are available,
though poor ones are still in common use. You may be happy with the
system-supplied random number generator on your computer, but you should
be aware that as computers get faster, requirements on random number
generators increase. Nowadays, a simulation that calls a random number
generator millions of times can often finish before you can make it down
the hall to the coffee machine and back.
A very nice review of random number generators was written by Pierre
L'Ecuyer, as Chapter 4 of the book: Handbook on Simulation, Jerry Banks,
ed. (Wiley, 1997). The chapter is available in postscript from from
L'Ecuyer's ftp site (see references). Knuth's volume on Seminumerical
Algorithms (originally published in 1968) devotes 170 pages to random
number generators, and has recently been updated in its 3rd edition
(1997). It is brilliant, a classic. If you don't own it, you should
stop reading right now, run to the nearest bookstore, and buy it.
A good random number generator will satisfy both theoretical and
statistical properties. Theoretical properties are often hard to obtain
(they require real math!), but one prefers a random number generator
with a long period, low serial correlation, and a tendency _not_ to
"fall mainly on the planes." Statistical tests are performed with
numerical simulations. Generally, a random number generator is used to
estimate some quantity for which the theory of probability provides an
exact answer. Comparison to this exact answer provides a measure of
"randomness".
File: gsl-ref.info, Node: The Random Number Generator Interface, Next: Random number generator initialization, Prev: General comments on random numbers, Up: Random Number Generation
The Random Number Generator Interface
=====================================
It is important to remember that a random number generator is not a
"real" function like sine or cosine. Unlike real functions, successive
calls to a random number generator yield different return values. Of
course that is just what you want for a random number generator, but to
achieve this effect, the generator must keep track of some kind of
"state" variable. Sometimes this state is just an integer (sometimes
just the value of the previously generated random number), but often it
is more complicated than that and may involve a whole array of numbers,
possibly with some indices thrown in. To use the random number
generators, you do not need to know the details of what comprises the
state, and besides that varies from algorithm to algorithm.
The random number generator library uses two special structs,
`gsl_rng_type' which holds static information about each type of
generator and `gsl_rng' which describes an instance of a generator
created from a given `gsl_rng_type'.
The functions described in this section are declared in the header
file `gsl_rng.h'.
File: gsl-ref.info, Node: Random number generator initialization, Next: Sampling from a random number generator, Prev: The Random Number Generator Interface, Up: Random Number Generation
The details of each generator are given later in this chapter.
- Random: void gsl_rng_set (const gsl_rng * R, unsigned long int S)
This function initializes (or `seeds') the random number
generator. If the generator is seeded with the same value of S on
two different runs, the same stream of random numbers will be
generated by successive calls to the routines below. If different
values of S are supplied, then the generated streams of random
numbers should be completely different. If the seed S is zero
then the standard seed from the original implementation is used
instead. For example, the original Fortran source code for the
`ranlux' generator used a seed of 314159265, and so choosing S
equal to zero reproduces this when using `gsl_rng_ranlux'.
- Random: void gsl_rng_free (gsl_rng * R)
This function frees all the memory associated with the generator R.
File: gsl-ref.info, Node: Sampling from a random number generator, Next: Auxiliary random number generator functions, Prev: Random number generator initialization, Up: Random Number Generation
Sampling from a random number generator
=======================================
The following functions return uniformly distributed random numbers,
either as integers or double precision floating point numbers. To
obtain non-uniform distributions *note Random Number Distributions::..
- Random: unsigned long int gsl_rng_get (const gsl_rng * R)
This function returns a random integer from the generator R. The
minimum and maximum values depend on the algorithm used, but all
integers in the range [MIN,MAX] are equally likely. The values of
MIN and MAX can determined using the auxiliary functions
This function returns a positive double precision floating point
number uniformly distributed in the range (0,1), excluding both
0.0 and 1.0. The number is obtained by sampling the generator
with the algorithm of `gsl_rng_uniform' until a non-zero value is
obtained. You can use this function if you need to avoid a
singularity at 0.0.
- Random: unsigned long int gsl_rng_uniform_int (const gsl_rng * R,
unsigned long int N)
This function returns a random integer from 0 to N-1 inclusive.
All integers in the range [0,N-1] are equally likely, regardless
of the generator used. An offset correction is applied so that
zero is always returned with the correct probability, for any
minimum value of the underlying generator.
If N is larger than the range of the generator then the function
calls the error handler with an error code of `GSL_EINVAL' and
returns zero.
File: gsl-ref.info, Node: Auxiliary random number generator functions, Next: Random number environment variables, Prev: Sampling from a random number generator, Up: Random Number Generation
Auxiliary random number generator functions
===========================================
The following functions provide information about an existing
generator. You should use them in preference to hard-coding the
These function return a pointer to the state of generator R and
its size. You can use this information to access the state
directly. For example, the following code will write the state of
a generator to a stream,
void * state = gsl_rng_state (r);
size_t n = gsl_rng_size (r);
fwrite (state, n, 1, stream);
File: gsl-ref.info, Node: Random number environment variables, Next: Saving and restoring random number generator state, Prev: Auxiliary random number generator functions, Up: Random Number Generation
Random number environment variables
===================================
The library allows you to choose a default generator and seed from
the environment variables `GSL_RNG_TYPE' and `GSL_RNG_SEED' and the
function `gsl_rng_env_setup'. This makes it easy try out different
generators and seeds without having to recompile your program.
File: gsl-ref.info, Node: Saving and restoring random number generator state, Next: Random number generator algorithms, Prev: Random number environment variables, Up: Random Number Generation
Saving and restoring random number generator state
This function prints a hex-dump of the state of the generator R to
`stdout'. At the moment its only use is for debugging.
File: gsl-ref.info, Node: Random number generator algorithms, Next: Unix random number generators, Prev: Saving and restoring random number generator state, Up: Random Number Generation
Random number generator algorithms
==================================
The functions described above make no reference to the actual
algorithm used. This is deliberate so that you can switch algorithms
without having to change any of your application source code. The
library provides a large number of generators of different types,
including simulation quality generators, generators provided for
compatibility with other libraries and historical generators from the
past.
The following generators are recommended for use in simulation. They
have extremely long periods, low correlation and pass most statistical
tests.
- Generator: gsl_rng_mt19937
The MT19937 generator of Makoto Matsumoto and Takuji Nishimura is a
variant of the twisted generalized feedback shift-register
algorithm, and is known as the "Mersenne Twister" generator. It
has a Mersenne prime period of 2^19937 - 1 (about 10^6000) and is
equi-distributed in 623 dimensions. It has passed the DIEHARD
statistical tests. It uses 624 words of state per generator and is
comparable in speed to the other generators. The original
generator used a default seed of 4357 and choosing S equal to zero
in `gsl_rng_set' reproduces this.
For more information see,
Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A
623-dimensionally equidistributed uniform pseudorandom number
generator". `ACM Transactions on Modeling and Computer
computed modulo 2^32. In the formulas above ^^ denotes
"exclusive-or". Note that the algorithm relies on the properties
of 32-bit unsigned integers and has been implemented using a
bitmask of `0xFFFFFFFF' to make it work on 64 bit machines.
The period of this generator is 2^88 (about 10^26). It uses 3
words of state per generator. For more information see,
P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe
Generators", `Mathematics of Computation', 65, 213 (1996),
203-213.
- Generator: gsl_rng_gfsr4
The `gfsr4' generator is like a lagged-fibonacci generator, and
produces each number as an `xor''d sum of four previous values.
r_n = r_{n-A} ^^ r_{n-B} ^^ r_{n-C} ^^ r_{n-D}
Ziff (ref below) notes that "it is now widely known" that two-tap
registers (such as R250, which is described below) have serious
flaws, the most obvious one being the three-point correlation that
comes from the defn of the generator. Nice mathematical
properties can be derived for GFSR's, and numerics bears out the
claim that 4-tap GFSR's with appropriately chosen offsets are as
random as can be measured, using the author's test.
This implementation uses the values suggested the the example on
p392 of Ziff's article: A=471, B=1586, C=6988, D=9689.
If the offsets are appropriately chosen (such the one ones in this
implementation), then the sequence is said to be maximal. I'm not
sure what that means, but I would guess that means all states are
part of the same cycle, which would mean that the period for this
generator is astronomical; it is (2^K)^D \approx 10^{93334} where
K=32 is the number of bits in the word, and D is the longest lag.
This would also mean that any one random number could easily be
zero; ie 0 <= r < 2^32.
Ziff doesn't say so, but it seems to me that the bits are
completely independent here, so one could use this as an efficient
bit generator; each number supplying 32 random bits.
For more information see,
Robert M. Ziff, "Four-tap shift-register-sequence
random-number generators", `Computers in Physics', 12(4),
Jul/Aug 1998, pp 385-392.
File: gsl-ref.info, Node: Unix random number generators, Next: Numerical Recipes generators, Prev: Random number generator algorithms, Up: Random Number Generation
Unix random number generators
=============================
The standard Unix random number generators `rand', `random' and
`rand48' are provided as part of GSL. Although these generators are
widely available individually often they aren't all available on the
same platform. This makes it difficult to write portable code using
them and so we have included the complete set of Unix generators in GSL
for convenience. Note that these generators don't produce high-quality
randomness and aren't suitable for work requiring accurate statistics.
However, if you won't be measuring statistical quantities and just want
to introduce some variation into your program then these generators are
quite acceptable.
- Generator: gsl_rng_rand
This is the BSD `rand()' generator. Its sequence is
x_{n+1} = (a x_n + c) mod m
with a = 1103515245, c = 12345 and m = 2^31. The seed specifies
the initial value, x_1. The period of this generator is 2^31, and
it uses 1 word of storage per generator.
- Generator: gsl_rng_random_bsd
- Generator: gsl_rng_random_libc5
- Generator: gsl_rng_random_glibc2
These generators implement the `random()' family of functions, a
set of linear feedback shift register generators originally used
in BSD Unix. There are several versions of `random()' in use
today: the original BSD version (e.g. on SunOS4), a libc5 version
(common on existing GNU/Linux systems) and a glibc2 version. Each
version uses a different seeding procedure, and thus produces
different sequences.
The original BSD routines accepted a variable length buffer for the
generator state, with longer buffers providing higher-quality
randomness. The `random()' function implemented algorithms for
buffer lengths of 8, 32, 64, 128 and 256 bytes, and the algorithm
with the largest length that would fit into the user-supplied
buffer was used. To support these algorithms additional
generators are available with the following names,
gsl_rng_random8_bsd
gsl_rng_random32_bsd
gsl_rng_random64_bsd
gsl_rng_random128_bsd
gsl_rng_random256_bsd
where the numeric suffix indicates the buffer length. The
original BSD `random' function used a 128-byte default buffer and
so `gsl_rng_random_bsd' has been made equivalent to
`gsl_rng_random128_bsd'. Corresponding versions of the `libc5'
and `glibc2' generators are also available, with the names
`gsl_rng_random8_libc5', `gsl_rng_random8_glibc2', etc.
- Generator: gsl_rng_rand48
This is the Unix `rand48' generator. Its sequence is
x_{n+1} = (a x_n + c) mod m
defined on 48-bit unsigned integers with a = 25214903917, c = 11
and m = 2^48. The seed specifies the upper 32 bits of the initial
value, x_1, with the lower 16 bits set to `0x330E'. The function
`gsl_rng_get' returns the upper 32 bits from each term of the
sequence. This does not have a direct parallel in the original
`rand48' functions, but forcing the result to type `long int'
reproduces the output of `mrand48'. The function
`gsl_rng_uniform' uses the full 48 bits of internal state to return
the double precision number x_n/m, which is equivalent to the
function `drand48'. Note that some versions of the GNU C Library
contained a bug in `mrand48' function which caused it to produce
different results (only the lower 16-bits of the return value were
set).
File: gsl-ref.info, Node: Numerical Recipes generators, Next: Other random number generators, Prev: Unix random number generators, Up: Random Number Generation
Numerical Recipes generators
============================
The following generators are provided for compatibility with
`Numerical Recipes'. Note that the original Numerical Recipes
functions used single precision while we use double precision. This
will lead to minor discrepancies, but only at the level of
single-precision rounding error. If necessary you can force the
returned values to single precision by storing them in a `volatile
float', which prevents the value being held in a register with double
or extended precision. Apart from this difference the underlying
algorithms for the integer part of the generators are the same.
- Generator: gsl_rng_ran0
Numerical recipes `ran0' implements Park and Miller's MINSTD
algorithm with a modified seeding procedure.
- Generator: gsl_rng_ran1
Numerical recipes `ran1' implements Park and Miller's MINSTD
algorithm with a 32-element Bayes-Durham shuffle box.
- Generator: gsl_rng_ran2
Numerical recipes `ran2' implements a L'Ecuyer combined recursive
generator with a 32-element Bayes-Durham shuffle-box.
File: gsl-ref.info, Node: Other random number generators, Next: Random Number Generator Performance, Prev: Numerical Recipes generators, Up: Random Number Generation
Other random number generators
==============================
The generators in this section are provided for compatibility with
existing libraries. If you are converting an existing program to use
GSL then you can select these generators to check your new
implementation against the original one, using the same random number
generator. After verifying that your new program reproduces the
original results you can then switch to a higher-quality generator.
Note that most of the generators in this section are based on single
linear congruence relations, which are the least sophisticated type of
generator. In particular, linear congruences have poor properties when
used with a non-prime modulus, as several of these routines do (e.g.
with a power of two modulus, 2^31 or 2^32). This leads to periodicity
in the least significant bits of each number, with only the higher bits
having any randomness. Thus if you want to produce a random bitstream
it is best to avoid using the least significant bits.
The following generator is provided for compatibility with the CRAY
MATHLIB routine RANF. It produces double precision floating point
numbers which should be identical to those from the original RANF.
- Generator: gsl_rng_ranf
This is the CRAY random number generator `RANF'. Its sequence is
x_{n+1} = (a x_n) mod m
defined on 48-bit unsigned integers with a = 44485709377909 and m
= 2^48. The seed specifies the lower 32 bits of the initial value,
x_1, with the lowest bit set to prevent the seed taking an even
value. The upper 16 bits of x_1 are set to 0. A consequence of
this procedure is that the pairs of seeds 2 and 3, 4 and 5, etc
produce the same sequences.
There is a subtlety in the implementation of the seeding. The
initial state is reversed through one step, by multiplying by the
modular inverse of a mod m. This is done for compatibility with
the original CRAY implementation.
Note that you can only seed the generator with integers up to
2^32, while the original CRAY implementation uses non-portable
wide integers which can cover all 2^48 states of the generator.
The function `gsl_rng_get' returns the upper 32 bits from each term
of the sequence. The function `gsl_rng_uniform' uses the full 48
bits to return the double precision number x_n/m.
The period of this generator is 2^46.
- Generator: gsl_rng_ranmar
This is the RANMAR lagged-fibonacci generator of Marsaglia, Zaman
and Tsang. It is a 24-bit generator, originally designed for
single-precision IEEE floating point numbers. It was included in
the CERNLIB high-energy physics library.
- Generator: gsl_rng_r250
This is the shift-register generator of Kirkpatrick and Stoll. The
sequence is
x_n = x_{n-103} ^^ x_{n-250}
where ^^ denote "exclusive-or", defined on 32-bit words. The
period of this generator is about 2^250 and it uses 250 words of
state per generator.
For more information see,
S. Kirkpatrick and E. Stoll, "A very fast shift-register
sequence random number generator", `Journal of Computational
Physics', 40, 517-526 (1981)
- Generator: gsl_rng_tt800
This is an earlier version of the twisted generalized feedback
shift-register generator, and has been superseded by the
development of MT19937. However, it is still an acceptable
generator in its own right. It has a period of 2^800 and uses 33
words of storage per generator.
For more information see,
Makoto Matsumoto and Yoshiharu Kurita, "Twisted GFSR
Generators II", `ACM Transactions on Modelling and Computer
Simulation', Vol. 4, No. 3, 1994, pages 254-266.
- Generator: gsl_rng_vax
This is the VAX generator `MTH$RANDOM'. Its sequence is,
x_{n+1} = (a x_n + c) mod m
with a = 69069, c = 1 and m = 2^32. The seed specifies the
initial value, x_1. The period of this generator is 2^32 and it
uses 1 word of storage per generator.
- Generator: gsl_rng_transputer
This is the random number generator from the INMOS Transputer
Development system. Its sequence is,
x_{n+1} = (a x_n) mod m
with a = 1664525 and m = 2^32. The seed specifies the initial
value, x_1.
- Generator: gsl_rng_randu
This is the IBM `RANDU' generator. Its sequence is
x_{n+1} = (a x_n) mod m
with a = 65539 and m = 2^31. The seed specifies the initial value,
x_1. The period of this generator was only 2^29. It has become a
textbook example of a poor generator.
- Generator: gsl_rng_minstd
This is Park and Miller's "minimal standard" MINSTD generator, a
simple linear congruence which takes care to avoid the major
pitfalls of such algorithms. Its sequence is,
x_{n+1} = (a x_n) mod m
with a = 16807 and m = 2^31 - 1 = 2147483647. The seed specifies
the initial value, x_1. The period of this generator is about
2^31.
This generator is used in the IMSL Library (subroutine RNUN) and in
MATLAB (the RAND function). It is also sometimes known by the
acronym "GGL" (I'm not sure what that stands for).
For more information see,
Park and Miller, "Random Number Generators: Good ones are
hard to find", `Communications of the ACM', October 1988,
Volume 31, No 10, pages 1192-1201.
- Generator: gsl_rng_uni
- Generator: gsl_rng_uni32
This is a reimplementation of the 16-bit SLATEC random number
generator RUNIF. A generalisation of the generator to 32 bits is
provided by `gsl_rng_uni32'. The original source code is
available from NETLIB.
- Generator: gsl_rng_slatec
This is the SLATEC random number generator RAND. It is ancient.
The original source code is available from NETLIB.
- Generator: gsl_rng_zuf
This is the ZUFALL lagged Fibonacci series generator of Peterson.
Its sequence is,
t = u_{n-273} + u_{n-607}
u_n = t - floor(t)
The original source code is available from NETLIB. For more
information see,
W. Petersen, "Lagged Fibonacci Random Number Generators for
the NEC SX-3", `International Journal of High Speed
Computing' (1994).
File: gsl-ref.info, Node: Random Number Generator Performance, Next: Random Number References and Further Reading, Prev: Other random number generators, Up: Random Number Generation
Random Number Generator Performance
===================================
The following table shows the relative performance of a selection the
available random number generators. The simulation quality generators
which offer the best performance are `taus', `gfsr4' and `mt19937'.
1754 k ints/sec, 870 k doubles/sec, taus
1613 k ints/sec, 855 k doubles/sec, gfsr4
1370 k ints/sec, 769 k doubles/sec, mt19937
565 k ints/sec, 571 k doubles/sec, ranlxs0
400 k ints/sec, 405 k doubles/sec, ranlxs1
490 k ints/sec, 389 k doubles/sec, mrg
407 k ints/sec, 297 k doubles/sec, ranlux
243 k ints/sec, 254 k doubles/sec, ranlxd1
251 k ints/sec, 253 k doubles/sec, ranlxs2
238 k ints/sec, 215 k doubles/sec, cmrg
247 k ints/sec, 198 k doubles/sec, ranlux389
141 k ints/sec, 140 k doubles/sec, ranlxd2
1852 k ints/sec, 935 k doubles/sec, ran3
813 k ints/sec, 575 k doubles/sec, ran0
787 k ints/sec, 476 k doubles/sec, ran1
379 k ints/sec, 292 k doubles/sec, ran2
File: gsl-ref.info, Node: Random Number References and Further Reading, Next: Random Number Acknowledgements, Prev: Random Number Generator Performance, Up: Random Number Generation
References and Further Reading
==============================
The subject of random number generation and testing is reviewed
extensively in Knuth's `Seminumerical Algorithms'.
Donald E. Knuth, `The Art of Computer Programming: Seminumerical
Algorithms' (Vol 2, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896842.
Further information is available in the review paper written by Pierre
